## 1. LOAD PACKAGES, DATA -----
rm(list=ls())
setwd(here::here())
# setwd() #Note: set your working directory here.
getwd()

library(tidyverse)
library(stringr)
library(sf)
library(lubridate)
library(ggExtra)
library(estimatr)
library(sensemakr)
library(interflex)
library(ggeffects)
library(ggpubr)
library(texreg)
library(sensemakr)
library(xtable)
library(cobalt)
library(MatchIt)
library(lmtest)
library(sandwich)
import::from(magrittr, "%<>%")
import::from(magrittr, "extract2")

df <- read.csv("hunnicutt_replication_analysis_data.csv", stringsAsFactors = F, strip.white = T)

df_matched <- matchit(pko_dep.d ~ predeployment_conflict + unmil_instrument + dist_adm2.cap + 
                        road_density_2000 + pop_density2000 + ntl_2003 + dist_gold + esa_forest,
                      data = df %>% 
                        select(ID_3, date, pko_dep.d, 
                               predeployment_conflict, road_density_2000, pop_density2000, 
                               unmil_instrument, dist_adm2.cap, ntl_2003, dist_gold, esa_forest) %>% 
                        filter_at(vars(predeployment_conflict, road_density_2000, pop_density2000, 
                                       unmil_instrument, dist_adm2.cap, ntl_2003, dist_gold, esa_forest), ~!is.na(.)),
                      method = "cem") %>% 
  match.data() %>% 
  select(ID_3, date, weights) %>% 
  left_join(df)

df_concessions <- st_read("./hunnicutt_replication_concessions_data/hunnicutt_replication_concessions_data.shp")

df_unmil_bases <- read.csv("hunnicutt_replication_unmil_bases_data.csv") %>% 
  st_as_sf(coords=c("longitude", "latitude"),
           crs=4326)

## 2. GENERATE TABLES -----
## ** Table 1: UN police are positively associated with the onset of natural resource concessions. -----
m1 <- lm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct +
           predeployment_conflict + unmil_instrument + 
           dist_adm2.cap +
           adj.c_new + adj.l1.c_new +
           road_density_2000 + pop_density2000 +
           ntl_2003 + 
           dist_gold + esa_forest +
           p1 + p2 + l1.conflict + conflict_roll.mean + 
           wb_projects + 
           qtr_yr + 
           t1 + t2 + t3, 
         data=df_matched, weights=weights) #main specification, ols
se1 <- coeftest(m1, vcov=NeweyWest(m1, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
p1 <- coeftest(m1, vcov=NeweyWest(m1, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

m2 <- glm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
            predeployment_conflict + unmil_instrument + 
            dist_adm2.cap +
            adj.c_new + adj.l1.c_new +
            road_density_2000 + pop_density2000 +
            ntl_2003 + 
            dist_gold + esa_forest +
            p1 + p2 + l1.conflict + conflict_roll.mean + 
            wb_projects + 
            t1 + t2 + t3, 
          data=df_matched, weights=weights, family="binomial") #main specification, logit
se2 <- coeftest(m2, vcov=NeweyWest(m2, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
p2 <- coeftest(m2, vcov=NeweyWest(m2, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

m3 <- lm_robust(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
                  predeployment_conflict + unmil_instrument + #covariates
                  dist_adm2.cap +
                  adj.c_new + adj.l1.c_new +
                  road_density_2000 + pop_density2000 +
                  ntl_2003 + 
                  dist_gold + esa_forest +
                  p1 + p2 + l1.conflict + conflict_roll.mean + 
                  wb_projects + 
                  qtr_yr + 
                  t1 + t2 + t3, 
                data=df_matched, weights=weights, clusters=as.factor(ID_1)) #main specification + adm1-clustered SEs, ols
se3 <- m3$std.error
p3 <- m3$p.value

m4 <- lm_robust(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
                  predeployment_conflict + unmil_instrument + 
                  dist_adm2.cap +
                  adj.c_new + adj.l1.c_new +
                  road_density_2000 + pop_density2000 +
                  ntl_2003 + 
                  dist_gold + esa_forest +
                  p1 + p2 + l1.conflict + conflict_roll.mean + 
                  wb_projects + 
                  qtr_yr + 
                  t1 + t2 + t3, 
                data=df_matched, weights=weights, clusters=as.factor(ID_2)) #main specification + adm2-clustered SEs, ols
se4 <- m4$std.error
p4 <- m4$p.value

m5 <- lm(c_new.d ~ l1.unpol_dep.ct + l12.unpol_dep.ct +
           l1.untrp_dep.ct  + l12.untrp_dep.ct +
           predeployment_conflict + unmil_instrument + 
           dist_adm2.cap +
           adj.c_new + adj.l1.c_new +
           road_density_2000 + pop_density2000 +
           ntl_2003 + 
           dist_gold + esa_forest +
           p1 + p2 + l1.conflict + conflict_roll.mean + 
           wb_projects + 
           qtr_yr + 
           t1 + t2 + t3, 
         data=df_matched, weights=weights) #main specification + 1-yr. lags of troops and police, ols
se5 <- coeftest(m5, vcov=NeweyWest(m5, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
p5 <- coeftest(m5, vcov=NeweyWest(m5, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

m6 <- lm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
           predeployment_conflict + unmil_instrument + 
           dist_adm2.cap +
           adj.c_new + adj.l1.c_new +
           road_density_2000 + pop_density2000 +
           ntl_2003 + 
           dist_gold + esa_forest +
           p1 + p2 + l1.conflict + conflict_roll.mean + 
           wb_projects + 
           qtr_yr + 
           t1 + t2 + t3, 
         data=df_matched %>% filter(l1.untrp_dep.ct<=0.15), weights=weights) #main specification + low-troop subgroup, ols
se6 <- coeftest(m6, vcov=NeweyWest(m6, lag = (4 * (nrow(df_matched %>% filter(l1.untrp_dep.ct<=0.15))/100)^(2/9)), prewhite = FALSE))[,2]
p6 <- coeftest(m6, vcov=NeweyWest(m6, lag = (4 * (nrow(df_matched %>% filter(l1.untrp_dep.ct<=0.15))/100)^(2/9)), prewhite = FALSE))[,4]

m7 <- lm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
           ldv_new_1mo + ldv_new_12mo + 
           predeployment_conflict + unmil_instrument + 
           dist_adm2.cap +
           adj.c_new + adj.l1.c_new +
           road_density_2000 + pop_density2000 +
           ntl_2003 + 
           dist_gold + esa_forest +
           p1 + p2 + l1.conflict + conflict_roll.mean + 
           wb_projects + 
           qtr_yr + 
           t1 + t2 + t3, 
         data=df_matched, weights=weights) #main specification + lagged dv, ols
se7 <- coeftest(m7, vcov=NeweyWest(m7, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
p7 <- coeftest(m7, vcov=NeweyWest(m7, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

m8 <- lm(c_new.d ~ l1.unpol_dep.ct * l1.untrp_dep.ct + 
           predeployment_conflict + unmil_instrument + 
           dist_adm2.cap +
           adj.c_new + adj.l1.c_new +
           road_density_2000 + pop_density2000 +
           ntl_2003 + 
           dist_gold + esa_forest +
           p1 + p2 + l1.conflict + conflict_roll.mean + 
           wb_projects + 
           qtr_yr + 
           t1 + t2 + t3, 
         data=df_matched, weights=weights) #police troops interaction, ols
se8 <- coeftest(m8, vcov=NeweyWest(m8, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
p8 <- coeftest(m8, vcov=NeweyWest(m8, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

texreg(list(m1, m3, m4, m5, m6, m8, m7, m2), 
       override.se = list(se1, se3, se4, se5, se6, se8, se7, se2), 
       override.pvalues = list(p1, p3, p4, p5, p6, p8, p7, p2),
       omit.coef = "t1|t2|t3|qtr_",
       include.ci=FALSE,
       include.rmse=FALSE,
       include.loglik=FALSE,
       include.deviance=FALSE,
       include.bic=FALSE,
       include.rsquared=FALSE,
       custom.coef.names = c("Intercept",
                             "UN Police$_{t-1}$ (100s)", "UN Troops$_{t-1}$ (1000s)",
                             "Pre-Deployment Conflict", "Distance to First UNMIL Bases",
                             "Distance to Nearest District Capital",
                             "Adjacent Investment", "Adjacent Investment$_{t-1}$",
                             "Pre-Deployment Road Density", "Pre-Deployment Population Density",
                             "Pre-Deployment Nighttime Luminosity", "Distance to Gold Deposit",
                             "Average Forest Cover", "Timber Sanctions Lifted", "Diamond Sanctions Lifted",
                             "Local Conflict$_{t-1}$", "Local Conflict (rolling mean)",
                             "World Bank Projects",
                             "UN Police$_{t-12}$ (100s)", "UN Troops$_{t-12}$ (1000s)",
                             "Police X Troops",
                             "New Investment$_{t-1}$", "New Investment$_{t-12}$"),
       reorder.coef = c(2, 19, 3, 20, 21, 18, 16, 17, 14, 15, 7, 8, 22, 23, 4, 5, 6, 9:13, 1),
       custom.header = list("Estimator: OLS"=1:7, "Estimator: GLM"=8),
       custom.model.names = paste0("(", letters[1:8], ")"),
       custom.gof.rows = list(`Time Trends`=rep("Yes", 8),
                              `Fixed Effects`=c(rep("Yes", 7), "No"),
                              `Clustered SEs`=c("", "County", "District", "", "", "", "", "")),
       caption = "UN police are positively associated with the onset of natural resource concessions.",
       label = "tab:results_summary")

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Table F1: Covariate Balance, Matching -----
matchit(pko_dep.d ~ predeployment_conflict + unmil_instrument + dist_adm2.cap + 
          road_density_2000 + pop_density2000 + ntl_2003 + dist_gold + esa_forest,
        data = df %>% 
          select(ID_3, date, pko_dep.d, 
                 predeployment_conflict, road_density_2000, pop_density2000, 
                 unmil_instrument, dist_adm2.cap, ntl_2003, dist_gold, esa_forest) %>% 
          filter_at(vars(predeployment_conflict, road_density_2000, pop_density2000, 
                         unmil_instrument, dist_adm2.cap, ntl_2003, dist_gold, esa_forest), ~!is.na(.)),
        method = "cem") %>% 
  bal.tab(un = T, stats = c("m", "v", "ks")) %>% 
  extract2("Balance") %>% 
  rownames_to_column(., var = "var") %>% 
  mutate(var = c("Pre-UNMIL Conflict", "Distance to 2004 Bases", "Distance to Nearest District Capital",
                 "Road Density", "Population Density", "Nighttime Lights", "Distance to Nearest Gold Deposit", "ESA Forest Cover")) %>% 
  select(var, Diff.Un, Diff.Adj, V.Ratio.Un, V.Ratio.Adj) %>% 
  xtable() %>% 
  print(include.rownames = F)

## ** Table F2: Additional Robustness Checks -----
rm1 <- lm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
            predeployment_conflict + unmil_instrument + 
            dist_adm2.cap +
            adj.c_new + adj.l1.c_new +
            road_density_2000 + pop_density2000 +
            ntl_2003 +
            dist_gold + esa_forest +
            p1 + p2 + l1.conflict + conflict_roll.mean +
            wb_projects +
            qtr_yr + 
            t1 + t2 + t3, 
          data=df_matched %>%
            filter(paste(ID_3, date) %in% c("180 2016-10-01", "80 2011-01-01", "84 2011-01-01", "84 2011-02-01") == F),
          weights=df_matched %>%
            filter(paste(ID_3, date) %in% c("180 2016-10-01", "80 2011-01-01", "84 2011-01-01", "84 2011-02-01") == F) %>%
            pull(weights)) #main specification + excluding outliers, ols
rm1_stats <- coeftest(rm1, vcov=NeweyWest(rm1, lag = 4*((7124/100)^(2/9)), prewhite = FALSE))
rse1 <- rm1_stats[,2]
rp1 <- rm1_stats[,4]

rm2 <- glm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
             predeployment_conflict + unmil_instrument + 
             dist_adm2.cap +
             adj.c_new + adj.l1.c_new +
             road_density_2000 + pop_density2000 +
             ntl_2003 + 
             dist_gold + esa_forest +
             p1 + p2 + l1.conflict + conflict_roll.mean + 
             wb_projects + 
             t1 + t2 + t3, 
           data=df_matched %>%
             filter(paste(ID_3, date) %in% c("180 2016-10-01", "80 2011-01-01", "84 2011-01-01", "84 2011-02-01") == F),
           weights=df_matched %>%
             filter(paste(ID_3, date) %in% c("180 2016-10-01", "80 2011-01-01", "84 2011-01-01", "84 2011-02-01") == F) %>%
             pull(weights), 
           family="binomial") #main specification + excluding outliers, glm
rm2_stats <- coeftest(rm2, vcov=NeweyWest(rm2, lag = 4*((7124/100)^(2/9)), prewhite = FALSE))
rse2 <- rm2_stats[,2]
rp2 <- rm2_stats[,4]

rm3 <- lm(c_new_intl.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
            predeployment_conflict + unmil_instrument + 
            dist_adm2.cap +
            adj.c_new + adj.l1.c_new +
            road_density_2000 + pop_density2000 +
            ntl_2003 + 
            dist_gold + esa_forest +
            p1 + p2 + l1.conflict + conflict_roll.mean + 
            wb_projects +
            qtr_yr + 
            t1 + t2 + t3, 
          data=df_matched %>% 
            mutate(c_new_intl.d = ifelse(c_new_intl > 0, 1, NA),
                   c_new_intl.d = ifelse(c_new_intl == 0, 0, c_new_intl.d)),
          weights=weights) #main specification + excluding unidentifiable agreements, ols
rse3 <- coeftest(rm3, vcov=NeweyWest(rm3, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
rp3 <- coeftest(rm3, vcov=NeweyWest(rm3, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

rm4 <- glm(c_new_intl.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
             predeployment_conflict + unmil_instrument + 
             dist_adm2.cap +
             adj.c_new + adj.l1.c_new +
             road_density_2000 + pop_density2000 +
             ntl_2003 + 
             dist_gold + esa_forest +
             p1 + p2 + l1.conflict + conflict_roll.mean + 
             wb_projects +
             t1 + t2 + t3, 
           data=df_matched %>% 
             mutate(c_new_intl.d = ifelse(c_new_intl > 0, 1, NA),
                    c_new_intl.d = ifelse(c_new_intl == 0, 0, c_new_intl.d)),
           weights=weights,
           family="binomial") #main specification + excluding unidentifiable agreements, glm
rse4 <- coeftest(rm4, vcov=NeweyWest(rm4, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
rp4 <- coeftest(rm4, vcov=NeweyWest(rm4, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

texreg(list(rm1, rm3, rm2, rm4), 
       override.se = list(rse1, rse3, rse2, rse4), 
       override.pvalues = list(rp1, rp3, rp2, rp4),
       omit.coef = "t1|t2|t3|qtr_",
       include.ci=FALSE,
       include.rmse=FALSE,
       include.loglik=FALSE,
       include.deviance=FALSE,
       include.bic=FALSE,
       include.rsquared=FALSE,
       custom.coef.names = c("Intercept",
                             "UN Police$_{t-1}$ (100s)", "UN Troops$_{t-1}$ (1000s)",
                             "Pre-Deployment Conflict", "Distance to First UNMIL Bases",
                             "Distance to Nearest District Capital",
                             "Adjacent Investment", "Adjacent Investment$_{t-1}$",
                             "Pre-Deployment Road Density", "Pre-Deployment Population Density",
                             "Pre-Deployment Nighttime Luminosity", "Distance to Gold Deposit",
                             "Average Forest Cover", "Timber Sanctions Lifted", "Diamond Sanctions Lifted",
                             "Local Conflict$_{t-1}$", "Local Conflict (rolling mean)",
                             "World Bank Projects"),
       reorder.coef = c(2:18, 1),
       custom.header = list("Estimator: OLS"=1:2, "Estimator: GLM"=3:4),
       custom.model.names = paste0("(", letters[1:4], ")"),
       custom.gof.rows = list(`Time Trends`=rep("Yes", 4),
                              `Fixed Effects`=c("Yes", "Yes", "No", "No"),
                              `Dropped Outliers`=c("Yes", "No", "Yes", "No"),
                              `Excluded Unidentifiable Concession Agreements`=c("No", "Yes", "No", "Yes")),
       caption = "Additional Robustness Checks",
       label = "si_tab:additional_robustness")

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Table F3: Omitted Variable Bias Sensitivity Analysis -----
sensemakr(lm(c_new.d ~ `UN Police (t-1)` + l1.untrp_dep.ct + 
               predeployment_conflict + unmil_instrument + 
               dist_adm2.cap +
               adj.c_new + adj.l1.c_new +
               road_density_2000 + pop_density2000 +
               ntl_2003 + 
               gold + esa_forest +
               p1 + p2 + l1.conflict + conflict_roll.mean + 
               wb_projects + 
               qtr_yr + 
               t1 + t2 + t3, 
             data=df_matched %>% rename(`UN Police (t-1)` = "l1.unpol_dep.ct", gold="dist_gold"), 
             weights=weights), 
          treatment="`UN Police (t-1)`",
          benchmark_covariates="gold",
          kd=c(3, 6, 9)) %>% 
  ovb_minimal_reporting()

## ** Table F4: UN police are not associated with the onset of riots, mob violence, or vigilante violence. -----
riots <- lm(riots.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
              predeployment_conflict + unmil_instrument + 
              dist_adm2.cap +
              adj.c_new + adj.l1.c_new +
              road_density_2000 + pop_density2000 +
              ntl_2003 + 
              dist_gold + esa_forest +
              p1 + p2 + l1.conflict + conflict_roll.mean + 
              qtr_yr + 
              t1 + t2 + t3, 
            data=df_matched, 
            weights=weights)
riots_se <- coeftest(riots, vcov=NeweyWest(riots, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
riots_p <- coeftest(riots, vcov=NeweyWest(riots, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

mob <- lm(mob.violence.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
            predeployment_conflict + unmil_instrument + 
            dist_adm2.cap +
            adj.c_new + adj.l1.c_new +
            road_density_2000 + pop_density2000 +
            ntl_2003 + 
            dist_gold + esa_forest +
            p1 + p2 + l1.conflict + conflict_roll.mean + 
            qtr_yr + 
            t1 + t2 + t3, 
          data=df_matched, weights=weights)
mob_se <- coeftest(mob, vcov=NeweyWest(mob, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
mob_p <- coeftest(mob, vcov=NeweyWest(mob, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

vigilante <- lm(vigilantes.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
                  predeployment_conflict + unmil_instrument + 
                  dist_adm2.cap +
                  adj.c_new + adj.l1.c_new +
                  road_density_2000 + pop_density2000 +
                  ntl_2003 + 
                  dist_gold + esa_forest +
                  p1 + p2 + l1.conflict + conflict_roll.mean + 
                  qtr_yr + 
                  t1 + t2 + t3,
                data=df_matched, weights=weights)
vigilante_se <- coeftest(vigilante, vcov=NeweyWest(vigilante, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
vigilante_p <- coeftest(vigilante, vcov=NeweyWest(vigilante, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

texreg(list(riots, mob, vigilante), 
       override.se = list(riots_se, mob_se, vigilante_se), 
       override.pvalues = list(riots_p, mob_p, vigilante_p),
       omit.coef = "t1|t2|t3|qtr_",
       include.ci=FALSE,
       include.rmse=FALSE,
       include.loglik=FALSE,
       include.deviance=FALSE,
       include.bic=FALSE,
       include.rsquared=FALSE,
       custom.coef.names = c("Intercept",
                             "UN Police$_{t-1}$ (100s)", "UN Troops$_{t-1}$ (1000s)",
                             "Pre-Deployment Conflict", "Distance to First UNMIL Bases", 
                             "Distance to Nearest District Capital",
                             "Adjacent Investment", "Adjacent Investment$_{t-1}$",
                             "Pre-Deployment Road Density", "Pre-Deployment Population Density",
                             "Pre-Deployment Nighttime Luminosity", "Distance to Gold Deposit",
                             "Average Forest Cover", "Timber Sanctions Lifted", "Diamond Sanctions Lifted",
                             "Local Conflict$_{t-1}$", "Local Conflict (rolling mean)"),
       reorder.coef = c(2:17, 1),
       custom.header = list("Outcome: Riots (0/1)"=1, "Outcome: Mob Violence (0/1)"=2, "Outcome: Vigilante Violence (0/1)"=3),
       custom.gof.rows = list(`Time Trends`=rep("Yes", 3),
                              `Fixed Effects`=c(rep("Yes", 3))),
       caption = "UN police are not associated with the onset of riots, mob violence, or vigilante violence.",
       label = "si_tab:mech1_resuls")

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Table F5: UN police are positively associated with short-term concession agreements but not long-term concession agreements. -----
short1 <- lm(dv ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
               predeployment_conflict + unmil_instrument + 
               dist_adm2.cap +
               adj.c_new + adj.l1.c_new +
               road_density_2000 + pop_density2000 +
               ntl_2003 + 
               dist_gold + esa_forest +
               p1 + p2 + l1.conflict + conflict_roll.mean + 
               wb_projects +
               qtr_yr + 
               t1 + t2 + t3, 
             data=df_matched %>% 
               mutate(dv = ifelse(c_new_short3yrs > 0, 1, NA),
                      dv = ifelse(c_new_short3yrs == 0, 0, dv)),
             weights=weights)
short1se <- coeftest(short1, vcov=NeweyWest(short1, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
short1p <- coeftest(short1, vcov=NeweyWest(short1, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

short2 <- lm(dv ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
               predeployment_conflict + unmil_instrument + 
               dist_adm2.cap +
               adj.c_new + adj.l1.c_new +
               road_density_2000 + pop_density2000 +
               ntl_2003 + 
               dist_gold + esa_forest +
               p1 + p2 + l1.conflict + conflict_roll.mean + 
               wb_projects +
               qtr_yr + 
               t1 + t2 + t3,
             data=df_matched %>% 
               mutate(dv = ifelse(c_new_short1q > 0, 1, NA),
                      dv = ifelse(c_new_short1q == 0, 0, dv)),
             weights=weights)
short2se <- coeftest(short2, vcov=NeweyWest(short2, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
short2p <- coeftest(short2, vcov=NeweyWest(short2, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

short3 <- lm(dv ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
               predeployment_conflict + unmil_instrument + 
               dist_adm2.cap +
               adj.c_new + adj.l1.c_new +
               road_density_2000 + pop_density2000 +
               ntl_2003 + 
               dist_gold + esa_forest +
               p1 + p2 + l1.conflict + conflict_roll.mean + 
               wb_projects +
               qtr_yr + 
               t1 + t2 + t3, 
             data=df_matched %>% 
               mutate(dv = ifelse(c_new_shortmean > 0, 1, NA),
                      dv = ifelse(c_new_shortmean == 0, 0, dv)),
             weights=weights)
short3se <- coeftest(short3, vcov=NeweyWest(short3, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
short3p <- coeftest(short3, vcov=NeweyWest(short3, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

long1 <- lm(dv ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
              predeployment_conflict + unmil_instrument + 
              dist_adm2.cap +
              adj.c_new + adj.l1.c_new +
              road_density_2000 + pop_density2000 +
              ntl_2003 + 
              dist_gold + esa_forest +
              p1 + p2 + l1.conflict + conflict_roll.mean + 
              wb_projects +
              qtr_yr + 
              t1 + t2 + t3,
            data=df_matched %>% 
              mutate(dv = ifelse(c_new_long10yrs > 0, 1, NA),
                     dv = ifelse(c_new_long10yrs == 0, 0, dv)),
            weights=weights)
long1se <- coeftest(long1, vcov=NeweyWest(long1, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
long1p <- coeftest(long1, vcov=NeweyWest(long1, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

long2 <- lm(dv ~ l1.unpol_dep.ct + l1.untrp_dep.ct +
              predeployment_conflict + unmil_instrument +
              dist_adm2.cap +
              adj.c_new + adj.l1.c_new +
              road_density_2000 + pop_density2000 +
              ntl_2003 + 
              dist_gold + esa_forest +
              p1 + p2 + l1.conflict + conflict_roll.mean + 
              wb_projects +
              qtr_yr +
              t1 + t2 + t3,
            data=df_matched %>% 
              mutate(dv = ifelse(c_new_long3q > 0, 1, NA),
                     dv = ifelse(c_new_long3q == 0, 0, dv)),
            weights=weights)
long2se <- coeftest(long2, vcov=NeweyWest(long2, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
long2p <- coeftest(long2, vcov=NeweyWest(long2, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

long3 <- lm(dv ~ l1.unpol_dep.ct + l1.untrp_dep.ct +
              predeployment_conflict + unmil_instrument +
              dist_adm2.cap +
              adj.c_new + adj.l1.c_new +
              road_density_2000 + pop_density2000 +
              ntl_2003 + 
              dist_gold + esa_forest +
              p1 + p2 + l1.conflict + conflict_roll.mean + 
              wb_projects +
              qtr_yr +
              t1 + t2 + t3,
            data=df_matched %>% 
              mutate(dv = ifelse(c_new_longmean > 0, 1, NA),
                     dv = ifelse(c_new_longmean == 0, 0, dv)),
            weights=weights)
long3se <- coeftest(long3, vcov=NeweyWest(long3, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,2]
long3p <- coeftest(long3, vcov=NeweyWest(long3, lag = (4 * (nrow(df_matched)/100)^(2/9)), prewhite = FALSE))[,4]

texreg(list(short1, short2, short3, long1, long2, long3), 
       override.se = list(short1se, short2se, short3se, long1se, long2se, long3se), 
       override.pvalues = list(short1p, short2p, short3p, long1p, long2p, long3p),
       omit.coef = "t1|t2|t3|qtr_",
       include.ci=FALSE,
       include.rmse=FALSE,
       include.loglik=FALSE,
       include.deviance=FALSE,
       include.bic=FALSE,
       include.rsquared=FALSE,
       custom.coef.names = c("Intercept",
                             "UN Police$_{t-1}$ (100s)", "UN Troops$_{t-1}$ (1000s)",
                             "Pre-Deployment Conflict", "Distance to First UNMIL Bases", 
                             "Distance to Nearest District Capital",
                             "Adjacent Investment", "Adjacent Investment$_{t-1}$",
                             "Pre-Deployment Road Density", "Pre-Deployment Population Density",
                             "Pre-Deployment Nighttime Luminosity", "Distance to Gold Deposit",
                             "Average Forest Cover", "Timber Sanctions Lifted", "Diamond Sanctions Lifted",
                             "Local Conflict$_{t-1}$", "Local Conflict (rolling mean)", "World Bank Projects"),
       reorder.coef = c(2:18, 1),
       custom.header = list("Short-term Concession Agreements"=1:3, "Long-term Concession Agreements"=4:6),
       custom.model.names = paste0("(", letters[1:6], ")"),
       custom.gof.rows = list(`Time Trends`=rep("Yes", 6),
                              `Fixed Effects`=c(rep("Yes", 6))),
       caption = "UN police are positively associated with short-term concession agreements but not long-term concession agreements.",
       label = "si_tab:contract_length")

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Table F6: High/Low Rule of Law Subgroup Results -----
med <- names(df_matched)[grepl("\\_a2", names(df_matched))]
spec <- c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct +
  predeployment_conflict + unmil_instrument + 
  dist_adm2.cap +
  adj.c_new + adj.l1.c_new +
  road_density_2000 + pop_density2000 +
  ntl_2003 + 
  dist_gold + esa_forest +
  p1 + p2 + l1.conflict + conflict_roll.mean + 
  wb_projects + 
  qtr_yr + 
  t1 + t2 + t3
sub_results <- list()
for(i in 1:length(med)){
  grp <- df_matched %>% #subgroup analysis
    group_by(ID_3) %>%
    summarise(med=mean(!!as.name(med[i]), na.rm=T), med=ifelse(is.nan(med), NA, med)) %>%
    mutate(med_qt = as.numeric(gtools::quantcut(med,q=2, na.rm=T))) %>% filter(!is.na(med_qt))  
  low_clan <- grp %>% filter(med_qt==1) %>% pull(ID_3)
  high_clan <- grp %>% filter(med_qt==2) %>% pull(ID_3)
  
  low <- lm(spec, data=df_matched %>% filter(ID_3%in%low_clan), weights=weights)
  low <- coeftest(low, 
                  vcov=NeweyWest(low, prewhite=FALSE, 
                                 lag = (4 * (nrow(df_matched[df_matched$ID_3%in%low_clan,])/100)^(2/9))))
  low <-  broom::tidy(low, conf.int=T) %>%
    filter(term%in%c("l1.untrp_dep.ct", "l1.unpol_dep.ct")) %>% 
    mutate(med=med[i], 
           grp="low",
           nobs=nobs(lm(spec, data=df_matched %>% filter(ID_3%in%low_clan), weights=weights)))
  
  high <- lm(spec, data=df_matched %>% filter(ID_3%in%high_clan), weights=weights)
  high <- coeftest(high, 
                   vcov=NeweyWest(high, prewhite=FALSE, 
                                  lag = (4 * (nrow(df_matched[df_matched$ID_3%in%high_clan,])/100)^(2/9))))
  high <-  broom::tidy(high, conf.int=T) %>%
    filter(term%in%c("l1.untrp_dep.ct", "l1.unpol_dep.ct")) %>% 
    mutate(med=med[i], 
           grp="high",
           nobs=nobs(lm(spec, data=df_matched %>% filter(ID_3%in%high_clan), weights=weights)))
  
  sub_results[[i]] <- rbind(low,high)
}
sub_results <- do.call("rbind", sub_results)

main_sub <- sub_results %>%
  filter(med%in%c("courts_eval_a2", "police_eval_a2", "fear_crime_a2", "red_crime_a2") & term=="l1.unpol_dep.ct") %>%
  mutate(se_tab=paste("(", round(std.error, digits=3), ")", sep=""),
         est_tab=ifelse(p.value>0.1, round(estimate, digits=3), NA),
         est_tab=ifelse(p.value<=0.1, paste(round(estimate, digits=3), "$^{+}$", sep=""), est_tab),
         est_tab=ifelse(p.value<=0.05, paste(round(estimate, digits=3), "$^{*}$", sep=""), est_tab),
         est_tab=ifelse(p.value<=0.01, paste(round(estimate, digits=3), "$^{**}$", sep=""), est_tab),
         est_tab=ifelse(p.value<=0.001, paste(round(estimate, digits=3), "$^{***}$", sep=""), est_tab),
         coef_lab=ifelse(term=="l1.unpol_dep.ct", "UN Police$_{t-1}$", "UN Troops$_{t-1}$")) #working here, 2.1
main_sub_t <- data.frame(var=c("UN Police$_{t-1}$", NA, "Fixed-Effects?", "Covariates?", "N"),
                         courts_eval_l=c(main_sub %>% filter(med=="courts_eval_a2" & grp=="low") %>% pull("est_tab"),
                                         main_sub %>% filter(med=="courts_eval_a2" & grp=="low") %>% pull("se_tab"),
                                         "Qtr-Yr", "Yes", main_sub %>% filter(med=="courts_eval_a2" & grp=="low") %>% pull("nobs")),
                         courts_eval_h=c(main_sub %>% filter(med=="courts_eval_a2" & grp=="high") %>% pull("est_tab"),
                                         main_sub %>% filter(med=="courts_eval_a2" & grp=="high") %>% pull("se_tab"),
                                         "Qtr-Yr", "Yes", main_sub %>% filter(med=="courts_eval_a2" & grp=="high") %>% pull("nobs")),
                         police_eval_l=c(main_sub %>% filter(med=="police_eval_a2" & grp=="low") %>% pull("est_tab"),
                                         main_sub %>% filter(med=="police_eval_a2" & grp=="low") %>% pull("se_tab"),
                                         "Qtr-Yr", "Yes", main_sub %>% filter(med=="police_eval_a2" & grp=="low") %>% pull("nobs")),
                         police_eval_h=c(main_sub %>% filter(med=="police_eval_a2" & grp=="high") %>% pull("est_tab"),
                                         main_sub %>% filter(med=="police_eval_a2" & grp=="high") %>% pull("se_tab"),
                                         "Qtr-Yr", "Yes", main_sub %>% filter(med=="police_eval_a2" & grp=="high") %>% pull("nobs")),
                         fear_crime_l=c(main_sub %>% filter(med=="fear_crime_a2" & grp=="low") %>% pull("est_tab"),
                                        main_sub %>% filter(med=="fear_crime_a2" & grp=="low") %>% pull("se_tab"),
                                        "Qtr-Yr", "Yes", main_sub %>% filter(med=="fear_crime_a2" & grp=="low") %>% pull("nobs")),
                         fear_crime_h=c(main_sub %>% filter(med=="fear_crime_a2" & grp=="high") %>% pull("est_tab"),
                                        main_sub %>% filter(med=="fear_crime_a2" & grp=="high") %>% pull("se_tab"),
                                        "Qtr-Yr", "Yes", main_sub %>% filter(med=="fear_crime_a2" & grp=="high") %>% pull("nobs")),
                         red_crime_l=c(main_sub %>% filter(med=="red_crime_a2" & grp=="low") %>% pull("est_tab"),
                                       main_sub %>% filter(med=="red_crime_a2" & grp=="low") %>% pull("se_tab"),
                                       "Qtr-Yr", "Yes", main_sub %>% filter(med=="red_crime_a2" & grp=="low") %>% pull("nobs")),
                         red_crime_h=c(main_sub %>% filter(med=="red_crime_a2" & grp=="high") %>% pull("est_tab"),
                                       main_sub %>% filter(med=="red_crime_a2" & grp=="high") %>% pull("se_tab"),
                                       "Qtr-Yr", "Yes", main_sub %>% filter(med=="red_crime_a2" & grp=="high") %>% pull("nobs")))
out <- list(main_sub_t)
attr(out, "message") <- c("Note: $^{+}$ p$<0.1$; $^{*}$ p$<0.05$; $^{**}$ p$<0.01$; $^{***}$ p$<0.001$")
print(xtableList(caption="High/Low Rule of Law Subgroup Results", out), include.rownames=F, sanitize.text.function=function(x){x})

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## 3. GENERATE FIGURES -----
## ** Figure 1a: Natural Resource Concessions in Post-Conflict Liberia - Sector -----
conc_ts_type <- data.frame(date=rep(seq.Date(from=as.Date("09-01-2004", "%m-%d-%Y"), to=as.Date("02-01-2018", "%m-%d-%Y"), by="month"), each=4),
                           count=NA, type=c("Agriculture", "Energy", "Forestry", "Mining"))
dates <- unique(conc_ts_type$date)
for(i in 1:length(dates)){
  d <- dates[i]
  conc_ts_type$count[conc_ts_type$date==d] <- c(sum(df_concessions$dt_strt<=d & df_concessions$sectors=="Agriculture"),
                                                sum(df_concessions$dt_strt<=d & df_concessions$sectors=="Energy Generation and Supply"),
                                                sum(df_concessions$dt_strt<=d & df_concessions$sectors=="Forestry"),
                                                sum(df_concessions$dt_strt<=d & df_concessions$sectors=="Industry, Mining, Construction"))
}

ggsave("./hunnicutt_replication_figures/figure1a.pdf", 
       ggplot(data=conc_ts_type) +
         geom_line(aes(date, count, col=type, group=type)) +
         theme_bw() +
         scale_x_date("Date") +
         scale_color_viridis_d(begin=0, end=0.75, guide=guide_legend(nrow=1)) +
         scale_y_continuous("Count of Concessions", limits=c(0, 360)) +
         theme(legend.position = "bottom", legend.title = element_blank(), text=element_text(size=8),
               legend.margin=margin(0,0,0,0),
               legend.box.margin=margin(-5,-5,-5,-5)), 
       width=6, height=3)

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Figure 1b: Natural Resource Concessions in Post-Conflict Liberia - Nationality of (Parent) Company -----
areas <- df_concessions %>% 
  st_drop_geometry() %>% 
  dplyr::select(aiddtID, sectors, company, country, sponsrs, pCompny, pContry) %>% 
  mutate(entity = NA,
         entity = ifelse(country == "Liberia" & !is.na(pContry), pContry, entity),
         entity = ifelse(country != "Liberia" & !is.na(country), country, entity),
         entity = ifelse(!is.na(country) & !is.na(pContry) & country != "Liberia", pContry, entity),
         entity = ifelse(is.na(entity) | country == "NA", "Missing information.", entity)) 

df_concessions %>% 
  st_drop_geometry() %>% 
  dplyr::select(aiddtID, sectors, company, country, sponsrs, pCompny, pContry) %>% 
  mutate(entity = NA,
         entity = ifelse(country == "Liberia" & !is.na(pContry), pContry, entity),
         entity = ifelse(country != "Liberia" & !is.na(country), country, entity),
         entity = ifelse(!is.na(country) & !is.na(pContry) & country != "Liberia", pContry, entity),
         entity = ifelse(is.na(entity) | country == "NA", "Missing information.", entity))  %>% 
  mutate(entity = ordered(factor(entity), levels = c("Missing information.", 
                                                     areas %>%
                                                       filter(grepl("information", entity)==F) %>%
                                                       group_by(entity) %>% 
                                                       summarise(total = n()) %>% 
                                                       arrange(total) %>%
                                                       pull(entity) %>% 
                                                       unique())),
         sectors = ifelse(sectors == "Energy Generation and Supply", "Energy", sectors),
         sectors = ifelse(sectors == "Industry, Mining, Construction", "Mining", sectors)) %>% 
  group_by(entity, sectors) %>% 
  summarise(total = n()) %>% 
  ggplot(., aes(entity, total, fill = sectors)) +
  geom_bar(stat = "identity", position = position_stack(reverse = TRUE)) +
  scale_fill_viridis_d(begin=0, end=0.75, guide=guide_legend(nrow=1)) +
  scale_y_continuous("Concessions Agreements Established") +
  scale_x_discrete("Nationality of Company or Parent Company") +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank(), text=element_text(size=8),
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(-5,-5,-5,-5))

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Figure 2a: UNMIL Bases near Guthrie Rubber Plantation, Bomi County, 2006 -----
bbox <- st_bbox(st_union(st_read("./LBR_adm_shp", "LBR_adm3") %>% 
                           st_transform(crs = 4326) %>% 
                           filter(NAME_1=="Bomi"))) %>% 
  st_as_sfc()
bomi_bases_2006 <- df_unmil_bases %>%
  filter(base %in% c("klay", "baha", "tubmanburg") & year==2006) %>%
  group_by(base) %>%
  summarise(pko=max(pko_deployed)) %>% 
  mutate(base = c("Baha", "Klay", "Tubmanburg"),
         pko = paste0("(PKO=", round(pko), ")"), 
         lab = paste(base, pko))
bomi_bases_2006 %<>%
  mutate(X = st_coordinates(bomi_bases_2006)[,1],
         Y = st_coordinates(bomi_bases_2006)[,2])

ggplot(data=NULL) +
  geom_sf(data=st_union(st_read("./LBR_adm_shp", "LBR_adm3") %>%
                          st_transform(crs = 4326)),
          fill="grey80", col="black", size=0.5) +
  geom_sf(data=df_concessions %>% filter(aiddtID==505), col="black", fill="grey30") +
  geom_sf(data=st_read("./grip4_liberia", "grip4_liberia") %>%
            st_set_crs(4236),
          col="grey50", size=0.05) +
  geom_sf(data=bomi_bases_2006, size=2) +
  geom_text(data=bomi_bases_2006 %>% filter(base=="Tubmanburg"), aes(X, Y, label = lab), colour = "black", hjust=-0.025, vjust=1.5, size=3.5) +
  geom_text(data=bomi_bases_2006 %>% filter(base=="Klay"), aes(X, Y, label = lab), colour = "black", hjust=-0.025, vjust=1.5, size=3.5) +
  geom_text(data=bomi_bases_2006 %>% filter(base=="Baha"), aes(X, Y, label = lab), colour = "black", hjust=-0.025, vjust=-1, size=3.5) +
  coord_sf(xlim = st_coordinates(bbox)[c(1,2),1], 
           ylim = st_coordinates(bbox)[c(2,3),2]) +
  ggspatial::annotation_scale() +
  theme_void()

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Figure 2b: UNMIL Bases near Cocopa Rubber Plantation, Nimba County, 2006 -----
bbox <- st_bbox(st_union(st_read("./LBR_adm_shp", "LBR_adm3") %>% 
                           st_transform(crs = 4326) %>% 
                           filter(NAME_2=="Saclepea"|NAME_3%in%c("Bain", "Garr")))) %>% 
  st_as_sfc()
nimba_bases_2006 <- df_unmil_bases %>%
  filter(base %in% c("ganta", "baila", "sagleipie") & year==2006) %>%
  group_by(base) %>%
  summarise(pko=max(pko_deployed)) %>% 
  mutate(base = c("Baila", "Ganta", "Sagleipie"),
         pko = paste0("(PKO=", round(pko), ")"), 
         lab = paste(base, pko))
nimba_bases_2006 %<>%
  mutate(X = st_coordinates(nimba_bases_2006)[,1],
         Y = st_coordinates(nimba_bases_2006)[,2])


ggplot(data=NULL) +
  geom_sf(data=st_union(st_read("./LBR_adm_shp", "LBR_adm3") %>%
                          st_transform(crs = 4326)),
          fill="grey80", col="black", size=0.5) +
  geom_sf(data=df_concessions %>% filter(aiddtID==451), col="black", fill="grey30") +
  geom_sf(data=st_read("./grip4_liberia", "grip4_liberia") %>%
            st_set_crs(4236),
          col="grey50", size=0.05) +
  geom_sf(data=nimba_bases_2006, size=2) +
  geom_text(data=nimba_bases_2006, aes(X, Y, label = lab), colour = "black", hjust=-0.025, vjust=1.5, size=3.5) +
  coord_sf(xlim = c(-9.192843, -8.65), 
           ylim = st_coordinates(bbox)[c(2,3),2]) +
  ggspatial::annotation_scale() +
  theme_void()

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])

## ** Figure 3a: Change in estimate given confounding X-times stronger -----
plot(sensemakr(lm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
                    predeployment_conflict + unmil_instrument +
                    dist_adm2.cap +
                    adj.c_new + adj.l1.c_new +
                    road_density_2000 + pop_density2000 +
                    ntl_2003 + 
                    gold + esa_forest +
                    p1 + p2 + l1.conflict + conflict_roll.mean + 
                    wb_projects + 
                    qtr_yr + 
                    t1 + t2 + t3,
                  data=df_matched %>% dplyr::rename(gold="dist_gold"), weights=weights),
               "l1.unpol_dep.ct", c("gold"), c(3, 6, 9)), 
     type="contour", sensitivity.of = "estimate", lim=0.3, lim.y=0.3)

## ** Figure 3b: Change in t-statistic given confounding X-times stronger -----
plot(sensemakr(lm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
                    predeployment_conflict + unmil_instrument +
                    dist_adm2.cap +
                    adj.c_new + adj.l1.c_new +
                    road_density_2000 + pop_density2000 +
                    ntl_2003 + 
                    gold + esa_forest +
                    p1 + p2 + l1.conflict + conflict_roll.mean + 
                    wb_projects + 
                    qtr_yr + 
                    t1 + t2 + t3,
                  data=df_matched %>% dplyr::rename(gold="dist_gold"), weights=weights),
               "l1.unpol_dep.ct", c("gold"), c(3, 6, 9)), 
     type="contour", sensitivity.of = "t-value", lim=0.3, lim.y=0.3)

## ** Figure A1: Distribution of Concession Areas per Agreement Type -----
data.frame(table(df_concessions$cntrctT)) %>%
  dplyr::rename(type="Var1", count="Freq") %>%
  mutate(type=str_replace_all(type, " \\s*\\([^\\)]+\\)", ""),
         type=fct_reorder(factor(type), count), 
         count_lab=paste0("n=", count)) %>%
  filter(count>0) %>%
  ggplot(.) +
  geom_bar(aes(type, count), stat="identity") +
  geom_text(aes(type, count, label=count_lab), hjust=-0.15, size=3.2) +
  scale_y_continuous(limits=c(0,200)) +
  labs(x="Type of Agreement", y="Count of Concession Areas") +
  theme_bw() +
  coord_flip()

## ** Figure E1: Largest Deployments of UNMIL Peacekeeping Police Outside of Monrovia -----
df_unmil_bases %>%
  filter(base %in% c("buchanan", "ganta", "gbarnga", "greenville", "harper", "kakata", "tapeta", "tubmanburg", "voinjama")) %>% 
  mutate(base=fct_reorder(factor(base), unpol, max),
         date=as.Date(date, "%Y-%m-%d"),
         base=str_to_title(base)) %>%
  ggplot(., aes(date, unpol, group=base)) +
  geom_line() +
  labs(x="Date", y="UNMIL Police Deployed") +
  theme_bw() +
  facet_wrap(.~base)


## ** Figure E2: Monthly Total of Peacekeepers Deployed to UNMIL, per clan -----
df %>% 
  mutate(pko_dep.ct=pko_dep.ct*1000, untrp_dep.ct=untrp_dep.ct*1000, unpol_dep.ct=unpol_dep.ct*100) %>% 
  filter(NAME_3 %in% unique(df$NAME_3[df$pko_dep.ct>0])) %>% 
  group_by(date, NAME_3) %>%
  summarise(pko=round(sum(pko_dep.ct)), untrp=round(sum(untrp_dep.ct)), unpol=round(sum(unpol_dep.ct))) %>%
  gather(., type, count, pko:unpol) %>%
  mutate(type=ordered(factor(type), levels=c("pko", "untrp", "unpol")),
         NAME_3=fct_reorder(NAME_3, count, min),
         date = as.Date(date, "%Y-%m-%d")) %>%
  ggplot(., aes(date, count, col=type)) +
  geom_line() +
  labs(x="Date", y="UNMIL Personnel Deployed") +
  scale_color_viridis_d("Personnel Type", labels=c("All Personnel", "Troops", "Police"), begin=0, end=0.75) +
  theme_bw() +
  theme(legend.position = "bottom", text=element_text(size=9)) +
  facet_wrap(.~NAME_3, ncol=11)

## ** Figure E3: Distribution of New Resource Concessions per Clan, 2004-2018 -----
ggplot(left_join(st_read("./LBR_adm_shp", "LBR_adm3") %>%
                   st_transform(crs = 4326), 
                 df %>% 
                   group_by(ID_3, year) %>%
                   summarise(new_conc=sum(c_new)))) +
  geom_sf(aes(fill=new_conc, group=ID_3), col="grey20", size=0.05) +
  scale_fill_gradient("New Concessions Established", low="white", high="black") +
  theme_void() +
  theme(legend.position = "bottom") +
  facet_wrap(.~year, ncol=5)

## ** Figure E4: Distribution of Active Resource Concessions per Clan, 2004-2018 -----
ggplot(left_join(st_read("./LBR_adm_shp", "LBR_adm3") %>%
                   st_transform(crs = 4326), 
                 df %>% 
                   group_by(ID_3, year) %>%
                   summarise(act_conc=sum(c_active)))) +
  geom_sf(aes(fill=act_conc, group=ID_3), col="grey20", size=0.05) +
  scale_fill_gradient("Active Concessions Established", low="white", high="black") +
  theme_void() +
  theme(legend.position = "bottom") +
  facet_wrap(.~year, ncol=5)

## ** Figure E7: Potential Outlying Clan-Months in Sample -----
main <- lm(c_new.d ~ l1.unpol_dep.ct + l1.untrp_dep.ct + 
             predeployment_conflict + unmil_instrument + 
             dist_adm2.cap +
             adj.c_new + adj.l1.c_new +
             road_density_2000 + pop_density2000 +
             ntl_2003 + 
             dist_gold + esa_forest +
             p1 + p2 + l1.conflict + conflict_roll.mean + 
             wb_projects +
             qtr_yr + 
             t1 + t2 + t3, 
           data=df_matched,
           weights=weights)

ggplot(data = NULL) +
  geom_point(data = df_matched %>% 
               filter_at(vars(c_new.d,
                              l1.unpol_dep.ct, l1.untrp_dep.ct, predeployment_conflict, unmil_instrument, 
                              dist_adm2.cap, adj.c_new, adj.l1.c_new, road_density_2000, pop_density2000,
                              ntl_2003, dist_gold, esa_forest, p1, p2, l1.conflict, conflict_roll.mean,
                              wb_projects, qtr_yr, t1, t2, t3), 
                         all_vars(!is.na(.))) %>% 
               ungroup() %>% 
               mutate(leverage_values = hatvalues(main),
                      residuals_norm = resid(main) / sqrt(sum(resid(main)^2)),
                      residuals_norm_sq= residuals_norm^2),
             aes(x = residuals_norm_sq, y = leverage_values),
             size = 1) + 
  geom_hline(yintercept = mean(hatvalues(main)),
             color = "red") +
  geom_vline(xintercept = mean(rstandard(main)), 
             color = "red") +
  ggrepel::geom_label_repel(data = df_matched %>% 
                              filter_at(vars(c_new.d,
                                             l1.unpol_dep.ct, l1.untrp_dep.ct, predeployment_conflict, unmil_instrument, 
                                             dist_adm2.cap, adj.c_new, adj.l1.c_new, road_density_2000, pop_density2000,
                                             ntl_2003, dist_gold, esa_forest, p1, p2, l1.conflict, conflict_roll.mean,
                                             wb_projects, qtr_yr, t1, t2, t3), 
                                        all_vars(!is.na(.))) %>% 
                              ungroup() %>% 
                              mutate(leverage_values = hatvalues(main),
                                     residuals_norm = resid(main) / sqrt(sum(resid(main)^2)),
                                     residuals_norm_sq= residuals_norm^2) %>%
                              filter(paste(ID_3, date) %in% c("180 2016-10-01", "80 2011-01-01", "84 2011-01-01", "84 2011-02-01")),
                            aes(x = residuals_norm_sq, y = leverage_values, label = paste0("Clan ID: ", ID_3, ",\n Date: ", date)),
                            min.segment.length = 0.2,
                            na.rm = TRUE,
                            size = 2) +
  labs(x = "Normalized residuals squared",
       y = "Leverage") +
  theme_minimal()

rm(list = ls()[ls() %in% c("df", "df_matched", "df_concessions", "df_unmil_bases") == F])






